package com.ruibang.glass.quality.util.spc;

/**
 * @Author: songJinKang
 * @CreateTime: 2023-10-18  15:50
 * @Description: TODO
 * @Version: 1.0
 */
public class ErrorFunction {
    private static final double EPSILON = 1e-6;
    private static final int MAX_ITERATIONS = 100;

    /**
     * @Description: 计算误差函数 erf(x) 的数值解
     * @version v1.0
     * @author songJinKang
     * @date 2023-10-18 16:02
     */
    public static double erf(double x) {
        double t = 1.0 / (1.0 + 0.5 * Math.abs(x));
        double ans = 1.0 - t *
                Math.exp(-x * x - 1.26551223 +
                        t * (1.00002368 +
                                t * (0.37409196 +
                                        t * (0.09678418 +
                                                t * (-0.18628806 +
                                                        t * (0.27886807 +
                                                                t * (-1.13520398 +
                                                                        t * (1.48851587 +
                                                                                t * (-0.82215223 +
                                                                                        t * 0.17087277)))))))));
        if (x < 0)
            return -ans;
        else {
            return ans;
        }
    }

    /**
     * @Description: 计算误差函数的反函数 erfinv(x) 的数值解
     * @version v1.0
     * @author songJinKang
     * @date 2023-10-18 16:01
     */
    public static double calculateInverseErf(double y) {
        if (y < -1.0 || y > 1.0) {
            return Double.NaN;
        } else if (y >= 0.0 && y < EPSILON) {
            return 0.0;
        } else if (y <= 0.0 && y > -EPSILON) {
            return 0.0;
        }

        double initGuess = Math.signum(y) * Math.sqrt(Math.sqrt(-Math.log(Math.min(1.0 - Math.abs(y), 1.0 - EPSILON))));
        double x = initGuess;

        for (int i = 0; i < MAX_ITERATIONS; i++) {
            double errorFunc = erf(x);
            double yPrime = Math.exp(-x * x) / Math.sqrt(Math.PI);
            double xNew = x + (y - errorFunc) / (yPrime * 2.0);
            if (x == xNew || Double.isNaN(xNew)) {
                break;
            }
            x = xNew;
        }

        return x;
    }

    /**
     * @Description: 计算 Cpk 对应的标准正态分布累积分布函数值
     * @version v1.0
     * @author songJinKang
     * @date 2023-10-18 16:01
     */
    public static double calculateZ(double cpk) {
        if (cpk < 0) {
            return Double.NaN;
        }
        double k = Math.log10(cpk * (Math.sqrt(1.0 - (cpk * cpk)) + 1.0));
        double a = 2.515517;
        double b = 0.802853;
        double c = 0.010328;
        double d = 1.432788;
        double e = 0.189269;
        double f = 0.001308;
        return Math.signum(cpk) * Math.sqrt(-2.0 * Math.log(Math.min(0.5 * Math.exp(-k * k) * (a + b * k + c * k * k + d * k * k * k + e * k * k * k * k + f * k * k * k * k * k), 1.0 - EPSILON)));
    }

}
